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We solve the rotational abrasion model of Roth, Marques and Durian [Phys. Rev. E (2010)], a 
one-dimensional quasilinear partial differential equation resembling the inviscid Burgers equation 
with the unusual feature of a step function factor as a coefficient. The complexity of the solution is 
primarily in keeping track of the cases in the piecewise function that results from certain amputation 
and interpolation processes, so we also extract from it a model of an evolving planar tree graph that 
tracks the evolution of the coarse features of the contour. 



What determines the shapes of pebbles is an intriguing 
physical question with interest not just to beachcombers 
out for walks but also geologists, who are interested in 
the history of erosion at a site [T] , as well as mechanical 
engineers [2], who wish to understand wear processes. 
Recently several models have been proposed to explain 
these shapes. Two stochastic models are of note, a "cut- 
ting model" O U| which accompanied an experimental 
measurement of pebbles rotating in a tray and an ana- 
lytically tractable "chipping model" [S]. These models 
lead to distributions of non-circular shapes. More re- 
cently, deterministic erosion processes have been stud- 
ied by Roth, Marques and Durian. They performed an 
experiment to measure the contours of linoleum tiles of 
fixed thickness and varying shape that they had rotated 
for differing amounts of time in a slurry of grit [6] . This 
paper describes the solution to their rotational abrasion 
model. Supposing r{9,t) describes the radial distance of 
the contour from the rotational axis as a function of angle 
6 and time t, they proposed 

Here C is a positive proportionality factor with dimen- 
sions of Angle/ (Time X Length) and H{x) denotes a Heav- 
iside step function defined so that H{x) = for a; < and 
H{x) = 1 for a; > 0. It is obvious that circles {r{6,t) = 
constant) are stationary solutions to this equation, and 
the evolution of experimental contours computed in [6] by 
a finite differencing scheme also evolved towards circles 
unerringly. These solutions also matched quantitatively 
the evolution of several geometrical quantities extracted 
from their experimental data, such as area, perimeter, 
and the width of the curvature distribution. 

Summarized here are the key ideas behind our exact 
solution of this equation. First, we exploit a connection 
to the Burgers equation at zero viscosity, a well-studied 
equation from gas dynamics [?]■ Second, the solution 
r{9,t) can be written in a piecewise fashion as a union 
of r = constant (circular) arcs and certain "stretched" 
segments of the initial contour. More precisely, these 
segments are curves of the form rQ{9{9Q,t)) where ro{9o) 
is the initial contour and 9{9o,t) at fixed 9q is a linear 
function in t. The solution is constructed to be con- 
tinuous, but will admit corners with discontinuous slope 
generically. Finally, the organization of the solution has 



a strong combinatorial flavor, and the evolution of the 
pattern of critical points in the contour is captured by a 
model of an evolving planar tree. This reduction suggests 
that discrete statistical models may capture the proper- 
tics of ensembles of abraded pebbles. The final section 
discusses other possible extensions. 

I. METHOD OF CHARACTERISTICS 

In what follows, we will usually think of the variables 
r and 9 in Eq. ([T]) as two dimensional rectangular coor- 
dinates (indeed, all the figures are plotted in this "un- 
wrapped" fashion) , though they do refer to polar coordi- 
nates, and we will usually refer to curves of constant r as 
circular arcs and to lengths in the 9 direction as angular 
widths. 

We first observe that Eq. ([T]) without the step func- 
tion factor is precisely the inviscid Burgers equation, a 
quasilinear first-order partial differential equation. That 
equation may be solved by the method of characteris- 
tics, which we shall now adapt. See also the book of 
Melikyan on solutions via characteristics to nonsmooth 
first-order equations in the theory of optimal control 
and in differential games [5- Let the initial contour 
be ro{9) = r{9,0) > 0. We now search for "charac- 
teristics", or space-time curves 9{a),t{a) beginning at 
9{Q) — 9q, t{0) — {a being some parameter) along 
which r{9(a),t{(7)) remains constant, and hence equal to 
ro(0o)- In other words, each point on the initial contour 
''o(^o) evolves forward in time along its characteristic. 

By applying the chain rule, we find 

±ri9icT),t{a))^0 

d9 dr dt dr ^ 
da 89 da dt 

Comparing this to Eq. ([T]), 



^ - CrH (^-^\ 
da~^'^[d9)- 

Integrating these equations with the initial condition 
t = 0, 9 = 9q and using the fact that r{9{a),t{a)) is 
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FIG. 1. (Color online) Evolution of a square, r vs9 (C=l) via 
the method of characteristics, t = (blue) and t = 0.04 (red). 
Dotted portions of the curve indicate pieces of the evolved 
contour which are amputated, and the horizontal dashes in- 
dicate discontinuities of the curve before amputation. The 
inset shows the two contours plotted in polar coordinates. 



points 9 > 9min ahead of this minimum will have traveled 
a nonzero angular distance along their characteristics at 
any t > 0, which results in a growing gap of undefined 
points (i.e. points such that Eq. ^ has no solution in 
Oq at a given t) between Omin and Omin + Crmint- The 
obvious thing to do is to interpolate by setting r = ?'mm 
for all in this interval, as this is the only natural way 
to ensure that the shape remains continuous. Thus inter- 
vals of constant r (circular arcs) are continually growing 
at local minima. 

The two cases just described are the simplest cases 
where the evolution along characteristics must be re- 
paired to become continuous. There are several more 
similar cases involving intervals of constant r which lead 
to multivaluedness (or no-valuedness) , but they are all 
treated by either amputation or interpolation, as de- 
scribed above. In the terminology of Melikyan [8], the 
points of amputation are "equivocal" and the points of 
interpolation are "dispersal" . 



constant and equal to ro{9a), we obtain 

rid{eo,t),t) = roieo) (2) 

e{eo,t) = eo + Cro{eo)tH(^^y (3) 

This equation yields a (possibly multi-valued) formal 
solution for r{9,t) via r{9{9o,t),t). To plot the evolu- 
tion of a curve with the solution in this form, begin with 
points distributed on the initial curve and move each of 
those points along its characteristic arc an angular dis- 
tance d9 — Cro{9o)H {d0r)dt in each time step dt. See 
Fig. [T] showing this evolution in the case of an initially 
square contour. 

II. MULTIVALUEDNESS AND AMPUTATION 

The solution given in the previous section is not yet 
well-defined; we must deal with the multivaluedness of 
the evolution along characteristics. To see how this 
arises, consider Fig. [T] Points on the contour slightly 
behind the local maxima will quickly overtake the points 
with the same r-values but slightly ahead of the maxi- 
mum, as those have nonpositive dgr and hence are frozen 
by the step function factor. This causes the two pieces of 
the contour to overlap, and is analogous to the formation 
of shocks in the inviscid Burgers equation. The "horizon- 
tal" discontinuity this creates is depicted in the figure as 
dotted lines. Based on the physical interpretation of the 
equation, the way to deal with this multivaluedness is to 
amputate the portion of the contour where this has oc- 
curred, as in the figure. This generates a corner in r{9, t) 
(though in the pictured example, the contour began with 
corners at the maxima). 

There are also problematic points around a local min- 
imum of r. Let irmim^min) be the coordinates of the 
local minimum on the initial contour. Since rmin > 0, all 



III. PIECEWISE SOLUTION FOR r{e,t) 

From the considerations above, giving an explicit for- 
mula for the solution would involve several layers of if- 
then constructs. We describe the full piecewise solution 
r{9, t) to Eq. ([T]) instead by decomposing the contour into 
strictly monotonic intervals, and within each of these the 
solution depends continuously on the initial contour. We 
also give the positions of the endpoints separating these 
intervals as a function of time. 

We define rising and falling faces of the contour to 
be connected components of points on the contour with 
dgr > and dgr < 0, respectively. See Fig. [2] for illustra- 
tions. The step function in Eq. ^ forces falling faces to 
be pieces of the original contour, i.e. r{9,t) = ro{9). The 
rising faces will be intervals from the original contour 
"sheared" by the evolution along characteristics. More 
precisely, r{9,t) is defined implicitly by solving for the 
value of ^0 in Eq. ^ such that 9 — 9{9Q,t), and then 
setting r{9, t) = ?'o(^. The shearing is caused by points 
moving with speed proportional to their radius. 

We will call intervals of constant r (circular) arcs^ and 
we will classify these into four types. Rising and falling 
arcs are those that are adjacent to rising, respectively 
falling faces on both sides. Min and max arcs are those 
which contain local minima, respectively maxima of the 
contour. Thus, r{9,t) at any fixed time decomposes into 
a set of faces and arcs. The evolution of each face or arc 
can be treated independently of the others for almost all 
times except for a discrete set of events when a face or 
arc changes into another type or disappears. 

The endpoints of faces or arcs fall into three categories 
named according to their behavior under time evolution: 
stationary endpoints, interpolating endpoints, and ampu- 
tating endpoints. Stationary endpoints are those that do 
not move under time evolution. There are two types, 
those at the right of a falling face and at the left of a 




falling arc 




FIG. 2. (Color online) Decomposition of the contour into faces 
and arcs; see text. The different types of endpoints are labeled 
with symbols: "0" are stationary, "+" are interpolating, "— " 
are amputating. 



min or falling arc, and those at the right of a max or 
falling arc and at the left of a falling face. Interpolating 
endpoints are those that move to the right under time 
evolution and are the sites of new interpolation. These 
are always to the right of rising or min arcs and to the left 
of rising faces. These endpoints move at constant speed 
Cr where r is the radius of the arc. Finally, amputating 
endpoints move to the right and are the sites of new am- 
putation. They are always to the right of rising faces and 
are on the left of max arcs, rising arcs, or falling faces. 

To calculate the motion of the amputating endpoints 
we will need the survival time of each point of the pebble. 
The intersection of the area under the contour with a cir- 
cle of radius r will be several disjoint arcs. From Eq. ([2| 
these shrink in time as the left endpoint of each moves 
with speed Cr towards the right, until the moment this 
value of r is amputated and the entire interval has van- 
ished. The remaining lifetime of a point (r, 9) inside the 



contour at time t is thus 



Cr 



where A{r,9,t) is the 




FIG. 3. (Color online) The face and arc merging processes, 
from top to bottom: vanishing of a max arc, vanishing of 
a falling face, and vanishing of a rising face. Endpoints are 
labeled as in Fig. [2] 



after the falling face determines which of the two faces 
vanish. When rr < rf, the falling face vanishes, when 
r,. > the rising face vanishes, and when — rf, they 
vanish simultaneously and the two arcs are joined. 

To make the above solution a bit more concrete, we 
sketch what happens in the case of a rectangular con- 
tour with side lengths 2a and 26 with a < b. Let 
7 = arctan(6/a). The initial contour (for Q < < 2tt) 
takes the form 



angular distance along the circle of radius r from (r, 0) ''o(^) 
to the left endpoint of its interval on the contour at time 
t. At fixed r, A(r, e,t) = O-Ol- Crt, where Ol is the 
position of this left endpoint on the initial contour, i.e. 
9l is an appropriate solution of ro{0L) — r. 

For an amputating endpoint that lies between a rising 
and falling face, the angular position as a function of time 
follows the contour of the falling face. Its position is a so- 
lution {r{t),e{t)) oiA{r,e,t) = Crt. One limit to keep in 
mind is when the falling face is vertical. Then the angu- 
lar position of the amputating endpoint will be stationary 
for a period of time. The opposite limit is when the am- 
putating endpoint lies between a rising face and an arc 
with radius r, then the endpoint moves to the right with 
constant speed Cr. The nature of the above decomposi- 
tion of the contour into faces and arcs changes precisely 
when endpoints collide with each other and faces and arcs 
merge. There are three basic processes: vanishing of a 
max arc, vanishing of a falling face, and vanishing of a 
rising face. See Fig. [3] for illustrations of these. If a rising 
face comes directly before a falling face, the radius of 
the arc before the rising face and the radius rf of the arc 



a sec < < 7 or 27r - 7 < < 27r 

b CSC 6 7<0<7r — 7 

-a sec 9 n — j<0<Tr + ^ 

-6 CSC 9 TT + 'y<9<27T — J. 



(4) 



In the first instant of time, the local minima at 9 — 
0, 7r/2, TT, 37r/2 expand by interpolation into min arcs, so 
the initial contour consists of a min arc of zero width 
at each of these local minima, each sandwiched between 
a falling face and a rising face. Fig. [l] depicts the first 
interval of time in the case b/a = 1, where the ampu- 
tating endpoints at the local maxima (initially at posi- 
tions 9 — 7, TT ± 7, 27r — 7) move to the right and the 
interpolating endpoints at the local minima do as well. 
The solution in this first time interval consists 12 piece- 
wise smooth curves, four each of min arcs, rising faces, 
and falling faces (though there is a tt periodicity of the 
contour in 9 which is preserved by the evolution, sim- 
plifying matters somewhat). This proceeds until time 
Ti — ^/^~'^^^^°'^(''/°) -y^fiien the points on the contour above 
r — b have all been amputated, and we remove the cor- 
responding rising and falling faces. If a ^ 6, the contour 
now consists of 8 piecewise smooth curves, two each of 
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FIG. 4. (Color online) (top) A schematic contour for a rect- 
angular contour after some abrasion and its associated tree. 
The two levels correspond to the distance Ri from the origin 
to the long sides and the distance R2 to the short sides. The 
lengths of the edges, corresponding to the widths of radial 
slices of the pebble, change in time according to their color 
and position - black edges at level i shrink with rate CRt, 
white edges sandwiched between black edges grow with rate 
CRi. (bottom) The contour and tree after the four original 
corners of the rectangle have been amputated. 

max arcs, min arcs, rising faces, and falling faces, as in 
the bottom left of Fig. ffl The max arcs vanish at time 

arcscc( — 6/a) — arcsccm/a) t-\ • n ^ ^ r 

T2 = — Durmg the hnal phase ot 

the evolution, the contour consists of 6 piecewise smooth 
curves, two each of min arcs, rising faces and falling 
faces. The rising faces and falling faces all vanish at time 
T3 — and for all later times the contour is a circle 
with radius r = a. 



IV. COARSE EVOLUTION AND TREE MODEL 

The constructions in the previous section are a bit un- 
wieldy to write out by hand, though it is straightforward 
to program a computer to map {0,t) to ro(0o) a-nd thus 
solve the evolution to the precision of the initial data 
ro{0o). The complexity is all in how the pattern of arcs 
and faces changes over time. In this section, we extract 
from the solution above a more combinatorial model of 
the evolution that focuses on this pattern and may make 
it more intuitive. 

First, identify all values of r such that the initial con- 
tour has a point where dgr — (do not count corners 
at maxima). Order these values from minimum to max- 
imum to define i?i < R2 < • ■ ■ < Rn- Geometrically, 
these critical values are the radii of circles centered at 
the origin which are tangent to the contour. These will 
be the N "levels" of a planar tree graph which we are 
constructing, which represents something like a skeleton 
of the contour. For example, regular polygons have only 
one level, the distance from the origin to any edge, and 
the rectangle discussed earlier has two levels, a and b. 

Consider the set of intervals arising from the intersec- 



tion of a circle of radius Ri with the area below the con- 
tour. Cut each interval in this set into black and white 
edges as follows: every subinterval which contains arcs 
(i.e. a segment on the contour coincident with the circle 
of radius Ri) becomes a white edge, and all other subin- 
tervals become black edges. Note that the black edges 
correspond to subintervals which support some hump of 
the contour. Each edge is assigned a length equal to 
the angular width of its subinterval. The edges just con- 
structed constitute all the edges at level i. An example 
of white and black edge assignments can be seen on the 
left of Fig. ^ 

We will now glue these edges into a planar tree. First, 
create a root vertex at level 1. Attach one end of every 
edge at level 1 to this root vertex, preserving the cyclic 
ordering of edges. Next, create a vertex at every black 
edge whose corresponding hump has intervals at level 2 
above it, and then attach one end of every edge at level 
2 to the appropriate vertex, preserving the linear order- 
ing. Repeat this process of creating vertices and gluing 
for each remaining level. See Fig. |4] for the example of a 
rectangle. Roughly speaking, the tree captures the pat- 
tern of the protrusions of the contour as we move from 
the origin outwards. Note that the edges at level i are 
not necessarily all attached to edges at level i — 1. 

The dynamics of endpoints, faces and arcs yields the 
following rules for the evolution of the tree. As time pro- 
gresses, every black edge at level i shrinks at the rate 
CRi. White edges at level i which happen to be sand- 
wiched between black edges (in the cyclic ordering around 
the root if i = 1, or the linear ordering above a vertex if 
« > 1) grow at the same rate CRi. White edges at level 
i > 1 which sit alone on a vertex shrink at the rate CRi. 
The lengths of all other white edges are held constant. 
If the length of an edge shrinks to zero, we remove it; if 
two white edges become adjacent on a vertex, we merge 
them into one white edge with length equal to the sum 
of their lengths. 

Under this evolution, the tree contracts from the leaves 
inwards; the edges supporting a branch will never vanish 
before the edges at higher levels connected to it. The to- 
tal time of evolution is thus determined by the length of 
the longest black edge at level 1. Translating back to the 
original contour, this means we just need to measure the 
angular width of the base of the largest hump. There- 
fore the contours with fixed minimum radius which take 
longest to evolve to a circle are those with a single min- 
imum radius. The tree picture also makes it clear that 
all contours evolve to a circle with radius equal to the 
minimum radius of the initial shape. 

At all times the pattern of white edges and black edges 
at different levels on the tree may be used to reconstruct a 
coarse version of the contour at that point in time. This is 
not a one-to-one correspondence between trees and con- 
tours, as there many possible contours that lead to the 
same tree. More explicitly, one can place at the left and 
right of the intervals corresponding to each black edge 
an arbitrary increasing function (rising face) between the 
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radii Ri to Ri+i, respectively decreasing function (falling 
face) , provided the angular widths of these two functions 
is consistent with the lengths of the edges above and the 
length of this black edge. Indeed, not even the maximum 
height of each hump enters this description. However, all 
contours that lead to the same tree have the same pat- 
tern of face and arc disappearances. This property gives 
some stability to the evolution - if noise is added to the 
initial contour, this will only affect the long term behav- 
ior insofar as it might change the widths of the bases 
of the large scale features. Small humps coming from 
the short-wavelength part of the noise will correspond to 
short edges on the tree which quickly vanish or merge 
with the large edges. 

V. DISCUSSION 

The solution in this paper generalizes easily to the case 
where the equation takes the form dtr + f{r)dgrH{der) = 
with nondecreasing / in place of Cr. If / is not nonde- 
creasing the solution above will be modified significantly 
as then some points of the contour would propagate in the 
opposite direction. However, a choice like this would also 
seem physically unmotivated. The function / allows for 
more general radius-speed relations, and the tree picture 
makes it clear that all that changes is the relative rate 
of growth or shrinkage of each edge, and not the overall 
qualitative picture; in particular, this may explain the 
observation in [6] that the model was robust to changing 
r to r" with a = 1/2, 1,2,3. 

We speculate next on some possible choices for /. The 
rectangle is the most interesting case studied by Roth, 
Marques and Durian, as it is the only contour with two 
widely-separated levels. In their data (see rightmost 
panel of Fig. 2 in [5]), two of the corners are abraded 
before the other two, whereas the model predicts that 
all four corners vanish at the same time (r2 in the nota- 



tion of the end of Sec. III). Thus the "constant" C may 
differ between the rising faces, perhaps being larger if 
the radii of the preceding minimum is smaller, and more 
generally, the speed / might in general be a function not 
just of r but also of r„iin of the face as well. In prin- 
ciple, (different branches of) /(r) can be extracted from 
experimentally measured curve contours by computing 
9tr, dgr and comparing them at fixed r, but preliminary 
analysis of data provided by Roth, Marques and Durian 
was not conclusive due to the difficulty of estimating dtr 
accurately. 

We did not carry out an extensive comparison here 
with the numerical solution of [S] , but the plotted curves 
appeared indistinguishable in a few checks. Indeed, it 
would be interesting to place this work on firmer math- 
ematical ground along the lines of ^8J by analyzing how 
the corners generated by the amputation and interpola- 
tion processes are smoothed by the addition of higher 
derivative terms, and how this happens in the finite dif- 
ference scheme of Roth, Marques and Durian. We also 
did not yet attempt to calculate the typical evolution 
of area, perimeter and other geometric quantities from 
our exact solution. Finally, we leave open the question 
whether erosion or shape evolution models in general may 
be simplified by posing them as laws for evolving tree 
graphs. In particular, it may be easier to construct mod- 
els for the evolution of an ensemble of pebbles in terms of 
a mean-field model on trees, rather than attempt a more 
direct description of ensembles of interacting contours. 
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